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ABSTRACT 

The crisis of the standard cooling flow model brought about by Chandra and XMM-Newton observa- 
tions of galaxy clusters, has led to the development of several models which explore different heating 
processes in order to assess if they can quench the cooling flow. Among the most appealing mechanisms 
are thermal conduction and heating through buoyant gas deposited in the ICM by AGNs. We combine 
Virgo/M87 observations of three satellites (Chandra, XMM-Newton and Beppo-SAX) to inspect the dy- 
namics of the ICM in the center of the cluster. Using the spectral deprojection technique, we derive the 
physical quantities describing the ICM and determine the extra-heating needed to balance the cooling 
flow assuming that thermal conduction operates at a fixed fraction of the Spitzer value. We assume that 
the extra-heating is due to buoyant gas and we fit the data using the model developed by Ruszkowski and 
Begelman (2002). We derive a scale radius for the model of ~ 5 kpc, which is comparable with the M87 
AGN jet extension, and a required luminosity of the AGN of a few x 10 42 ergs _1 , which is comparable to 
the observed AGN luminosity. We discuss a scenario where the buoyant bubbles are filled of relativistic 
particles and magnetic field responsible for the radio emission in M87. The AGN is supposed to be 
intermittent and to inject populations of buoyant bubbles through a succession of outbursts. We also 
study the X-ray cool component detected in the radio lobes and suggest that it is structured in blobs 
which are tied to the radio buoyant bubbles. 

Subject headings: conduction — cooling flows — galaxies: active — X-rays: galaxies: clusters — 
galaxies: clusters: individual (Virgo) 



1. INTRODUCTION 

The hot diffuse X-ray emitting gas (intracluster 
medium, ICM for short) provides a powerful tool to in- 
spect the internal dynamics of galaxy clusters. For the 
typical density and temperature of the intracluster gas, 
the main emission mechanism is the bremsstrahlung and, 
for a large amount of clusters, the radiative cooling time in 
the central regions is significantly shorter than the Hubble 
time. As a consequence, if no additional heating mech- 
anism is present, the gas cools and is expected to flow 
inwards, forming a cooling flow. The standard model of 
cooling flows (see Fabian 1994, for a review) predicted 
the gas to be a multiphase medium in which there is 
a broad range of temperatures and densities present at 
all radii. Mass deposition rates were estimated to be 
as large as hundreds of solar masses per year (Allen et 
al. 2001). This model was strengthened by the general 
thought that in presence of magnetic fields the thermal 
conduction must be highly suppressed (Binney and Cowie 
1981; Fabian 1994; Chandran and Cowley 1998; Malyshkin 
2001), which is a necessary condition for the multiphase 
cooling to operate. In fact, no heating exchange between 
the different phases must occur in order that they may 
coexist. There is some observational evidence that mod- 



est magnetic fields are present throughout the intracluster 
medium. The current measurements of intracluster mag- 
netic fields are based on Faraday rotation measure (RM) 
in radio sources seen through clusters (e.g. Kim, Kronberg 
and Tribble 1991; Clarke, Kronberg and Bohringer 2001; 
Feretti et al. 1999; Taylor et al. 2001); direct evidence 
also comes from measurements of extended regions of ra- 
dio synchrotron emission in clusters (see e.g. Giovannini 
and Feretti 2000; Fusco-Femiano et al. 2000; Owen, Mor- 
rison and Voges 1999; Feretti 1999). Both the excess RM 
values and the radio halo data suggest modest magnetic 
fields, at a few microgauss levels, throughout the cluster. 

Recent XMM-Newton and Chandra observations have 
shown that in the central regions, the temperature drops 
to about one third of its overall mean value and there is 
no evidence of temperatures smaller than ^1 — 2 keV 
(Peterson et al. 2001; Kaastra et al. 2001; Tamura et al. 
2001; Allen et al. 2001), suggesting that the gas does not 
cool below these cutoff temperatures. Moreover, the new 
estimated mass deposition rates are significantly smaller 
than those evaluated by using previous X-ray satellites 
data (McNamara et al. 2001; Peterson et al. 2001). Lastly, 
the new data show that clusters spectra are better repre- 
sented by a single (or double) temperature model rather 
than the standard multiphase (multi-temperature) cooling 
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flow model (Molendi and Pizzolato 2001; Bohringer et al. 
2001; Fabian et al. 2001; Matsushita et al. 2002). 

These new results clearly show that the standard cooling 
flow model is not a satisfying description of the internal 
dynamics of the ICM. Some source of heat which stops the 
cooling flow and balances radiative losses must be sought. 
The nature of this source and the origin of the heat mech- 
anism is still unclear. 

One possible candidate is thermal conduction. Recent 
works by Narayan and Medvedev (2001) and Gruzinov 
(2002) show that in the presence of turbulent magnetic 
fields, the conductivity can be as large as a fraction of the 
Spitzer value and thus can play a significant role in balanc- 
ing cooling flows. As a consequence, thermal conduction 
has been recently re- introduced as a possible heat source to 
balance the energy losses (see e.g. Voigt et al. 2002; Voigt 
and Fabian 2004; Fabian et al. 2002; Malyshkin 2001; Za- 
kamska and Narayan 2003). However, as we will discuss 
more in detail in §5.1, thermal conduction fails in supply- 
ing the needed heat in the central regions (see also Voigt et 
al. 2002; Voigt and Fabian 2004; Zakamska and Narayan 
2003). 

Heating from a central active galactic nucleus is an- 
other possibility. The idea is supported by the fact that 
most of "cooling clusters" host a central active galactic 
nucleus with strong radio activity (Burns 1990; Brighenti 
and Mathews 2002). Several models in the literature ex- 
plore AGN heating processes to assess if they can balance 
the radiative losses. One of the most appealing mecha- 
nisms involves buoyant gas bubbles, inflated by the AGN, 
that subsequently rise through the cluster ICM heating it 
up (Churazov et al. 2002; Bohringer et al. 2002; Briiggcn 
and Kaiser 2002a,b; Briiggen et al. 2002). 

However, all the models predicting that radiative cooling 
is balanced only by energy input from the central AGN fail 
(McNamara 2002; Zakamska and Narayan 2003; Brighenti 
and Mathews 2002). In particular, Brighenti and Math- 
ews (2002) analyzed several heating mechanisms induced 
by the central AGN and concluded that no simple mech- 
anism is able to quench the cooling flow. Moreover, the 
required mechanism needs a finely tuned heating source. 
Indeed, the heat source must provide sufficient energy to 
stop the cooling flow, but not enough to trigger strong con- 
vection or the metallicity gradients observed in all cooling 
flow clusters (De Grandi and Molendi 2001) would be de- 
stroyed. As a consequence, an AGN can be an efficient 
mechanism in the very center of the cluster but it is un- 
likely to be strong enough to provide energy to the outer 
parts of the "cooling region" . So, it may be viewed as com- 
plementary to thermal conduction which fails in quenching 
the cooling flow in the innermost regions. 

Recently, Ruszkowski and Begclman (2002) and Zakam- 
ska and Narayan (2003) concluded that both thermal con- 
duction and heating from a central AGN can play an 
important role in balancing the cooling. In particular, 
Ruszkowski and Begelman (RB02 hereafter) developed a 
model where both thermal conduction and heating from a 
central AGN co-operate in balancing the radiative losses. 
One of the main advantages of this model is that it reaches 
a stable final equilibrium state and it is able to reproduce 
the main observed quantities, such as the temperature pro- 
file (with a minimum temperature T ~ 1 keV) . 



In this paper, we use M87/Virgo observations of three 
satellites (namely XMM-Newton, Chandra and Beppo- 
SAX) to test various heating models on this cluster. To 
this end, we apply to the M87 data the deprojection tech- 
nique to recover some physical quantities of the ICM such 
as the gravitational mass, the entropy and the heating re- 
quired to balance the cooling flow. 

The paper is organized as follows: in §2 we report details 
about the analysis of the three (Chandra, XMM-Newton 
and Beppo-SAX) M87 datasets; in §3, we revise briefly 
the spectral deprojection technique that we adopt for our 
analysis; in §4, we deproject the M87 data, we test that 
the spherical symmetry hypothesis holds and we derive 
the gravitational mass for M87; in §5, we determine the 
amount of extra-heating needed to balance the cooling flow 
when thermal conductivity is assumed to operate at a frac- 
tion of the Spitzer value. Lastly, in §6 we summarize our 
results. 

2. DATA ANALYSIS 

Thanks to its proximity, the Virgo cluster and its giant 
elliptical central galaxy M87, represent an incomparable 
target to inspect the internal properties of the ICM. Aim- 
ing to a precise characterization of the ICM, we use ob- 
servations of the three satellites Chandra, XMM-Newton 
and Beppo-SAX. These satellites views are complemen- 
tary: the sharp PSF (~ 0.5") of Chandra can provide a 
precise analysis of the innermost (say ,$,5 — 10 kpc) re- 
gions of M87; the XMM-Newton large collecting area and 
its wide field of view allow a good inspection of the inter- 
mediate regions (up to <~ 80 kpc); in the outermost regions 
of the cluster, where the angular resolution is less critical 
and XMM-Newton data are highly contaminated by the 
background, Beppo-SAX is a better choice and data can 
be collected out to a radius of ~ 120 kpc. 

2.1. XMM-Newton data preparation. 

M87 has been observed during the PV phase of XMM- 
Newton. The details of this observation have been widely 
discussed in several publications (see e.g. Bohringer et al. 
2001; Belsole et al. 2001; Molendi and Pizzolato 2001; 
Molendi and Gastaldello 2001; Gastaldello and Molendi 
2002; Matsushita et al. 2002). We make use of the re- 
sults of the spectral analysis described and discussed in 
Molendi (2002, M02 hereafter). The cluster is divided in 
139 regions for 12 concentric annuli, centered on the emis- 
sion peak. The regions are the same as those presented 
in M02, apart from the annuli in the 1 — 4 arcmin range 
which have been taken 0.5 arcmin wide instead of 1 arcmin 
wide. Unlike Matsushita et al. (2002), we decided to use 
annuli at least 30" wide in order to avoid possible PSF 
contaminations (see also Markevitch 2002). For a detailed 
description of the MOS PSF see Ghizzardi (2001). 

We refer the reader to M02 for all the details of the 
spectral analysis procedure and remind the reader that 
the accumulated spectra in all the regions are fitted with 
two different models: (i) a single temperature (IT) model 
(vmekal in XSPEC) and (ii) a two temperature (2T) model 
(vmekal + vmekal in xspec). 

The IT fits of the accumulated spectra provide for each 
region the emission- weighted temperature T of the gas and 
the emission integral EI — J n e n p dV where n e and n p are 
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the electron and proton density. The 2T fits of the spectra 
provide for each region the temperatures and the emission 
integrals of the two different components of the gas. These 
quantities will be used to derive the density and the tem- 
perature profile of the cluster. 

2.2. Chandra data preparation. 

We have analyzed the Chandra ACIS-S3 observation 
(obs. id. 352; see also Young et al. 2002) centered on M87 
(a = 12:30:49.40; S = +12:23:27.82) using CIAO 2.1.1 and 
CALDB 2.15. We have followed the procedures described 
in the Science Threads available at the Chandra X-ray 
Center on-line pages. The light curve was filtered for high 
background events obtaining an effective exposure time of 
35.2 ks. From our analysis we excluded the AGN and the 
associated jet cutting off a narrow rectangular region cen- 
tered in a = 12:30:48.80, S = +12:23:31.37 (J2000) with 
lengths Aa = 25", A<5 = 5" and rotated by 22° from E to 
N. We have extracted spectra from regions of concentric 
annuli centered on the emission peak (0" — 10", 10" — 30", 
30"-45", 45"-60", 60"-80", 80"-100" and 100"-120"); 
the 10" — 30" annulus was divided into four 90° regions 
starting from a position angle of 45°, all the other annuli 
have been divided into 8 regions 45° wide starting from 
position angle 0°. The background used in the spectral 
fits was extracted from blank-sky observations using the 
acis-bkgrnd-lookup script. We have fitted each spectrum 
with the same models (IT and 2T) used for the analysis 
of the XMM-Newton spectra and using the effective areas 
and response matrix derived with the routines mkwarf and 
mkwrmf for extended sources given in CIAO. The energy 
range is the same as the one used for the XMM-Newton 
spectral analysis. 

2.3. Beppo-SAX data preparation. 

We have analyzed the pointed Beppo-SAX observation 
of the Virgo cluster (obs. id. 60010001), adding the 
MECS2 and MECS3 data and obtaining an effective ex- 
posure time of 25.1 ks. The analysis of data follows the 
procedure described in details in De Grandi and Molcndi 
(2001). We have extracted spectra for 7 concentric an- 
nuli centered on the emission peak (a — 187.6992 deg, 
5 = 12.3878 deg J2000), each annulus is 2' wide up to 10' 
from the peak and 4' wide from 10' to the maximum 20' 
radius. We have fitted each spectrum with a mekal model 
absorbed for the Galactic Nh using the appropriate ef- 
fective area computed for extended source as described in 
Dc Grandi and Molcndi (2001) and the background spec- 
trum extracted from blank-sky files for the same annular 
regions. The energy range considered in the spectral anal- 
ysis is 2. — 10. keV in all cases out of the 8' — 12' annulus for 
which we have considered the 3.5 — 10 keV range to avoid 
spectral distortions from the supporting structure of the 
instrument entrance windows. Note that, as we will dis- 
cuss in the next paragraph, for Beppo-SAX data, only the 
IT model has been used to fit the accumulated spectra. 

2.4. The joint data set. 

For each region observed with XMM-Newton or Chan- 
dra, an F-test is used to establish whether the 2T model 
provides a better description of the data with respect to 



the IT model. We used two different criteria for XMM- 
Newton and Chandra data to select the regions repre- 
sented by a 2T model. As far as XMM-Newton data are 
concerned, for those regions whose F-test provides a prob- 
ability > 95%, the 2T description is retained, whereas for 
those regions whose F-test provides a probability < 90%, 
the IT description is adopted. The regions having an F- 
test probability within the [0.90-0.95] range have been re- 
jected and excluded from the analysis. A somewhat more 
stringent criterion has been adopted for the Chandra data, 
regions with F-test probability within the [0.75-0.98] range 
have been rejected, because of the lower statistical quality 
of the latter dataset. 

It is worth noting that M02 shows (using XMM-Newton 
data) that the regions which are better represented by a 
2T model match the radio "arms" which are visible in the 
M87 map at 90 cm (Owen, Eilek and Kassim 2000). At 
large radii, where we have no evidence of a second tem- 
perature component from the XMM-Newton data, we will 
use the results of IT fits to the Beppo-SAX data, which, 
because of its limited spectral coverage, is insensitive to 
the cooler component. 

While the cool component is related to the radio "arms" , 
the hot component in the regions described by the 2T 
model, is very similar to the IT gas of the other regions 
which do not feature strong radio emission and are located 
at similar radial distances from the cluster core. In Fig. 
Al we plot the emission weighted temperatures for IT and 
2T regions. The open circles represent the temperatures 
of the IT regions, while the triangles plot the temperature 
of the hot component in those regions which are fitted by 
a 2T model. We plot error bars for a few representative 
points. All the other error bars are not reported for a clear 
reading of the plot. The Fig. Al shows that for those an- 
nuli where we have IT and 2T regions, the temperature 
of the hot component of the 2T regions is not separate 
from the temperature of the single phase gas of the IT re- 
gions. So, the hot component and the single phase gas are 
distributed in a regular and symmetric fashion. As far as 
only the hot component is considered for the 2T regions, 
the cluster appears to be approximativcly spherically sym- 
metric, which is an important condition for the application 
of the deprojection technique. 

As already outlined, for each region we consider the 
emission integral EI and the emission-weighted temper- 
ature T. For those regions fitted by the 2T model, we 
retain the EI and T related to the hot component. The 
radially averaged profile for each physical quantity is de- 
termined starting from the region-by-region description. 
We assign to each annulus, a mean T by averaging on all 
the T of the regions belonging to the annulus. The EI is 
the sum of the EI of the regions along the ring. 

Some attention must be paid in this averaging (or sum- 
mation for EI) procedure as some portions of the observed 
regions can be masked. In fact, there are some pixels of 
the Field of View of the different instruments which must 
be rejected for different reasons: (i) they are CCD hot 
or dead pixels; (ii) they correspond to the gaps between 
nearby CCDs; (iii) they correspond to regions which are 
never observed by the instruments (for MECS) or which 
are partially outside the Field of View; (iv) they corre- 
spond to some point source contaminating the X-ray clus- 
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ter emission; (v) they correspond to the excluded AGN 
rectangular regions {Chandra). As a consequence, even 
if the geometrical area (or equivalently the emitting vol- 
ume) of the regions belonging to the same ring is the same, 
the effective emitting area (or volume) is a fraction of the 
geometrical area depending on the number of rejected pix- 
els. The measured EI in each region accounts for photons 
coming from the effective emitting area (or volume); hence, 
each EI must be renormalized. The factor of normaliza- 
tion is given by the ratio A geom /(A geom - A mask ) where 
Ageom is the geometric area of the region and A mas k is the 
total area of the rejected pixels in the region. The same 
normalization factors are used as weights for the determi- 
nation of the averaged T. 

In Fig. A2, we plot the averaged profiles for T and 
El/Area, respectively in panels (a) and (b). Circles 
refer to XMM-Newton data, squares to Chandra data 
and triangles refer to Beppo-SAX data. The values for 
T and El/Area are also reported in Table Al. The 
emission integral is reported and plotted in xspec units, 
[1{)~ 1A / Aitd 2 ang (l + z) 2 ]EI where d ang is the angular dis- 
tance of the source in cm, z is the redshift and EI is in 
cm~ 3 . Area is the area of the ring in arcmin 2 . Note that 
error bars for the EI are quite small, especially for the 
XMM-Newton data thanks to its high effective area which 
allows a very precise measure of the emission integral up to 
^100 kpc. The EI j Area profiles from the three different 
data sets, match each other in the common ranges. On 
the contrary, Fig. A2(a) shows a systematic difference be- 
tween the three temperature profiles. The discrepancy be- 
tween XMM-Newton and Beppo-SAX is probably related 
to the use of different energy bands. For the observation 
of the Virgo cluster which has a temperature of 2.5 — 3 
keV, the XMM-Newton energy range (0.4 — 4 keV) is more 
suitable than the Beppo-SAX energy range (2 — 10 keV) 
and, consequently, XMM-Newton estimations are proba- 
bly more reliable. For what concerns the discrepancy be- 
tween XMM-Newton and Chandra temperature profiles, 
the differences are of the order of a few percent and well 
within the cross-calibration uncertainties between Chan- 
dra and XMM-Newton . To make full use of the combined 
Chandra, XMM-Newton and Beppo-SAX data and at the 
same time avoid unphysical jumps in deprojected quanti- 
ties when moving from one dataset to the next, we decided 
to shift, through a scale renormalization, the Chandra and 
the Beppo-SAX datasets. XMM-Newton was chosen as 
reference dataset because of its higher statistics. The scale 
factors have been derived by imposing that the three tem- 
perature profiles match in the common ranges. We find a 
renormalization factor x — 1.03 for the Chandra tempera- 
ture profile, and n = 1.10 for the Beppo-SAX temperature 
profile. 

The temperature profiles, corrected for the scale factor, 
are plotted in Fig. A3. The final joint data set used for 
the analysis is given by (i) the Chandra data in the 0' — 1' 
range (4 points); (ii) the XMM-Newton data in the 1' - 8' 
range (8 points); (iii) the Beppo-SAX data in the 8' — 12' 
range (3 points). 

3. SPECTRAL DEPROJECTION 

The deprojection technique has become very popular 
to investigate the intracluster medium properties (Ettori 



2002; Ettori et al. 2002; Matsushita et al. 2002; Allen ct 
al. 1996; Pizzolato et al. 2003). Under the assumption of 
spherical symmetry, the different 3-D quantities describing 
the ICM are derived from the 2-D projected ones, start- 
ing from the outer shell and moving inwards following an 
onion-peeling technique. 

Among the different prescriptions available for the de- 
projection, we decided to adopt the spectral deprojection 
introduced by Ettori et al. (2002) (see also Ettori, De 
Grandi and Molendi 2002). The physical quantities to 
be deprojected are those obtained from the spectral anal- 
ysis described in the previous Section. Each 3-D variable 
f shell is related to the projected one F ring according to the 
relation 

f shell = (V T ) #F ring , (1) 

where (V T ) 1 is the inverse of the transposed matrix V 
whose elements are the volumes of the i-th shell pro- 
jected on the j-th ring. The detailed evaluation of this ma- 
trix can be found in Kriss, Cioffi and Canizares (1983). By 
replacing F ring with (i) the emission-weighted measured 
temperature T r i ng L r i ng ; (ii) the ring luminosity L r i ng and 
(iii) (.E//0.82) 1 ' 2 , we can derive respectively eT s heii, the 
emissivity e and the electron density n e , where the relation 
EI = J n e n p dV = 0.82 J n 2 e dV has been used. The main 
advantage in using this technique is that deriving the elec- 
tron density n e from EI is straightforward without any 
assumption on its functional shape. Moreover, since our 
measurements of the EI are very accurate, an immediate 
and precise determination of the electron density profile 
can be obtained. It is worth noting that eq. (1) is derived 
using the onion-peeling procedure where the contribution 
to the emission in each shell is obtained from the projected 
quantity by subtracting off the emission contribution of 
the outer shells starting from the edge of the cluster and 
moving inwards. In addition to this basic prescription, a 
correction factor accounting for the cluster emission be- 
yond the maximum radius R max to infinity must be in- 
cluded. The procedure to evaluate this correction factor 
is presented in Appendix A. In practice, in eq. (1) the 
2-D variable to be deprojected (F ring ) is replaced by an 
effective one F^/J g (see eq. Al). 

From now on, in order to avoid confusion, we will use T 
for the 3-D deprojected temperature and Tew for the 2-D 
emission-weighted temperature. 

4. DEPROJECTING M87 

4.1. Deprojected profiles 

Following the prescription described in the previous Sec- 
tion, we derive e, n e and T profiles for M87. All the ex- 
tracted values are reported in Table A2. In Fig. A4 we 
plot (filled circles) the deprojected electron density n e and 
temperature T, respectively in panel (a) and (b). For com- 
parison, in Fig. A4 we overplot (open triangles) the depro- 
jected profiles obtained by Matsushita et al. (2002) who 
used a different choice of regions and a somewhat different 
technique for spectral deprojection. The Matsushita et al. 
(2002) results plotted here are those obtained by fitting 
the MOS data with a 2T model. Our profiles are in rea- 
sonable agreement with those derived by Matsushita et al. 
(2002). 
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Particular attention must be paid in the computation 
of the error bars for T and n e . In evaluating errors, we 
want to consider that, even if the deprojection technique 
assumes that the spherical symmetry condition is fulfilled, 
the dispersion of the EI and Tew measurements along 
each ring around the averaged value is often significantly 
larger than the error of each measure. In particular, this 
occurs for the XMM-Newton data, since the XMM-Newton 
large effective area allows a very good statistics for M87 
and provides EI measurements with very small error bars. 
The scatter of the data around the averaged value of the 
ring is a measure of the data displacement from the spher- 
ical symmetry. In order to account for this displacement 
in the final error evaluation, we decided to assign to the 
EI and Tew of each region, an error (<jei and ot bw re- 
spectively) which is the linear sum of two different contri- 
butions. The first contribution is simply the error derived 
from the spectral fit with the 1T/2T model. The second 
contribution is given by the dispersion of the measure- 
ments along the ring around the averaged value. Error 
bars for the deprojected quantities n e and T reported in 
Fig. A4 have been obtained running 1000 Monte-Carlo 
simulations, sampling the EI and Tew of each region 
around their mean value assuming Gaussian distributions 
for the errors <jei and <Jt E w ■ 

The profiles plotted in Fig. A4 will be used as start- 
ing points to derive some other quantities (such as mass 
and conductivity). In most cases, gradients of n e and 
T are involved. Some smoothing procedure will be re- 
quired to manage these derivatives. Consequently, we ap- 
ply a smoothing algorithm which replaces each point with 
the average value obtained using boxes of 3 points, i.e.: 
Vi = (Vi-i + Vi + Vj+i) /3. In smoothing the temperature 
profile, we excluded from the smoothing procedure the two 
last (Beppo-SAX) points, in order to preserve the final de- 
creasing behavior of the T profile. The final temperature 
profile is obtained by applying the smoothing procedure 
twice: firstly, we smooth the starting deprojected profile 
and then the obtained values are smoothed again. For 
the density profile, we applied the smoothing procedure 
separately to the first 6 points and the others in order 
to preserve the 2-/3 behavior. Again the two last points 
have been excluded from the smoothing operation and the 
smoothing has been applied twice. The open diamonds in 
Fig. A4 represent the n e and T profile after the smooth- 
ing operation. An alternative solution to smoothing is 
provided by the use of analytical functions fitting the pro- 
files. The solid lines in the figures are the best fits to the 
data. For the electron density profile we use the fitting 
function: 



n e (r) 



ni 
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(2) 



which corresponds to a 2-/3 model. For the temperature 
profile we find that the function 



T(r) =T -T 1 cxp - 



(3) 



provides a good description of the data. For the tem- 
perature profile, we find the following best-fit values: 
T = 2.399 ± 0.090 keV, T x = 0.776 ± 0.097 keV, a T = 



3.887' ±0.731'. For what concerns the electron density, we 
find n = 2.68' ± 0.54', a x = 0.71 ± 0.06, m = 0.033 ± 0.01 
cm" 3 , r 2 = 3.73' ± 9.65', a 2 = 20.19 ± 104.84 and 
n 2 = 0.069 ±0.010 cm -3 . The inferred best fit values have 
large statistical errors. This is due to the large number of 
free parameters adopted. So, we decided to fix two param- 
eters, namely the core radius and the slope of the second 
component; we fix the slope at large radii a 2 to the value 
obtained from the RASS measurements (Bohringer et al. 
1994) setting a 2 = 0.705 (/3 = 0.47). For the radius r 2 we 
decided to use the value of the ot inferred from the best fit 
of the temperature profile, which defines the scale radius 
for the rise of the temperature. Having fixed the values of 
r 2 and a 2 , we find: a x = 1.518 ±0.317, m = 0.089 ±0.011 
cm" 3 , n = 0.834' ± 0.175' and n 2 = 0.019 ± 0.002 cm" 3 . 
Note that the inferred value of n (<~ 5 kpc) roughly cor- 
responds to the AGN jet extension (e.g. Di Matteo et al. 
2003; Young et al. 2002). 

4.2. Sector deprojection 

The deprojection method is based on the assumption of 
spherical symmetry. As discussed in §2, the hot compo- 
nent in those regions which are described by the 2 tem- 
perature (2T) model behaves very similarly to the gas in 
the regions described by a single temperature (IT) model. 
However, we should like to verify if any correlation of the 
hot component with the radio emission exists, invalidat- 
ing our assumption of spherical symmetry. To this aim, we 
divided the cluster in sectors and deprojected separately 
each sector. Again we consider only the hot component for 
the 2T regions. The sectors have been chosen according to 
the radio emission regions: the [30° - 120°], [210° - 270°] 
and [330° — 360°] sectors have been cumulated together to 
form the non-radio sector. Furthermore, we analyzed sepa- 
rately each of the following sectors: [0°-30°], [120°-150°], 
[270° - 300°], [300° - 330°] which correspond or are close 
to the radio emission arms. 

The Chandra data are not suitable to perform such a 
study, because the statistics is not very high. In this 
case, possible azimuthal variations can be hidden by errors. 
Therefore, in order to verify the assumption of spherical 
symmetry, we consider only the XMM-Newton data (cir- 
cles in Fig. A2). Nevertheless, also with XMM-Newton 
data, it is quite difficult to use small sectors to perform 
sector-by-sector comparisons since their statistics is not 
very high. A significant comparison can be made between 
the whole cluster and the non-radio selected sector. In 
Fig. A5 (in panel (a) and (b) respectively) we compare 
the n e and T profiles for the whole cluster (filled dots) 
with the non-radio sector (open diamonds). In general, no 
significant differences arc evident. The profiles of the elec- 
tron density are almost identical whereas there is a slight 
tendency of the temperature calculated on the whole clus- 
ter to be smaller than the temperature of the non-radio 
sector. The difference is due to the fact that in the non- 
radio sector only IT regions contribute, while in the whole 
cluster profile also the contribution of the hot component 
temperature of the regions described by a 2T model is 
accounted. This temperature is slightly smaller than the 
overall temperature and produces a mild decrease of the 
whole cluster temperature profile with respect to the non- 
radio sector temperature profiles. However, differences arc 
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well within lcr; thus, excluding the radio emission sectors 
(or the regions described by the 2T model) does not affect 
significantly the results and no evident azimuthal asymme- 
try can be highlighted. We can conclude that the spherical 
symmetry assumption, which is an important condition for 
our analysis, is tenable. It is also worth noting that, in or- 
der to derive most of the other physical quantities (mass, 
conduction, entropy, etc.) a smoothing operation on T 
and n e is necessary, so that the small differences reported 
above would not anyhow affect our results. 

4.3. The gravitational mass for M87 

Once the basic quantities n e and T are obtained through 
the deprojection, under suitable assumptions, other re- 
lated quantities describing the ICM can be derived. One 
of the most important is the gravitational mass. Suppos- 
ing that the gas is in hydrostatic equilibrium within the 
potential well of the dark matter, the gravitational mass 
M within a radius r can be derived via the hydrostatic 
equilibrium equation: 



M(< r) = - 



kTr 
Gfirrip 



d In T d In ) 



din 5 



din'/ 



(4) 



where fx = 0.6 is the mean molecular weight, G the grav- 
itational constant and m p the proton mass. Even if eq. 
(4) provides a direct method to derive the gravitational 
mass, it is a differential equation and the temperature 
and the electron density are involved through their gra- 
dients. Irregular features in the profile induce jumps on 
the evaluated mass. A classic solution consists in smooth- 
ing the data and replacing the n e and the T profiles with 
their smoothed counterparts plotted in Fig. A4 (open dia- 
monds). As far as the temperature is concerned, a smooth- 
ing procedure is viable since errors are rather large and the 
general shape of the profile is quite smooth. Correspond- 
ingly, the smoothed profile is compatible (always within 
la) with the original one. On the contrary, for the electron 
density where errors are small, the smoothing procedure 
could hide some features which are physical. In order to 
assess whether the smoothing affects results, we compare 
the gravitational mass obtained using the temperature and 
the density smoothed profiles with the gravitational mass 
derived using the unsmoothed profiles where errors have 
been evaluated using the standard Monte Carlo technique. 
As we show in Fig. A6, the smoothed profile (filled cir- 
cles) agrees with the non-smoothed one (open diamonds) 
within 1 — 2(7 and no significant difference can be high- 
lighted. It is worth noting that the mass derived without 
smoothing provides three mass values (at r ~ 0.5, ~ 2 and 
<~ 25 kpc) which are negative and compatible with zero: 
M = -0.5t}i x 10 10 M Q , M = -2.7^ x 10 10 M Q and 
M = -Q.22±\^ x 10 12 M Q . For these points, in Fig. A6, 
we show only the upper limit of the error bar. 

Alternatively to the smoothing procedure, the analyti- 
cal expressions for n e and T (cqs. 2 and 3) can be used. 
The curve in Fig. A6 is the analytical mass obtained using 
these two best fit profiles. This mass profile has a plateaux 
at a radius of about 10—15 kpc (~ 2 arcmin). This flat- 
tening behavior is the consequence of the flattening of the 
2-/3 profile, at the same radius, which could correspond 
to the edge of the central cD. The analytical gravitational 



mass is very similar to the mass obtained both with the 
smoothing procedure and with the unsmoothed data. The 
differences with respect to the latter curve are limited to 
a few points (<~ 15, <~ 25 and ~ 90 kpc). For these points, 
some "holes" appear in the shape of the non-analytical 
profiles. It is worth noting that the "hole" in a (integrated) 
mass profile is not physical. It may indicate that in this 
region the hydrostatic equilibrium hypothesis breaks down 
(e.g. Pizzolato et al. 2003) and that an outflow providing 
additional pressure to support gravity occurs there (see 
next Section for equations and further details). In any 
case, the error bars of all these points are large enough to 
make "holes" compatible with the analytical profile and 
no strong evidence is present to claim that an outflow is 
present. 

For comparison, the grey-shaded regions in Fig. A6 
report the gravitational mass derived by Nulsen and 
Bohringcr (1995) using ROSAT-PSPC data. Our profile, 
in the common radial regions is in agreement. It is inter- 
esting to note that the point at r <~ 20 kpc has a very large 
error bar in both the estimations. 

5. COOLING, CONDUCTION AND HEATING 

As previously outlined, finding a heating mechanism 
able to balance the cooling is not an easy task. This 
Section will be devoted to the inspection of some heating 
sources using the M87 data set. 

The heating contribution required to balance the radia- 
tive cooling can be estimated starting from the thermody- 
namic equations describing a spherically symmetric clus- 
ter: 
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where <f> is the gravitational potential, M the gravitational 
mass within r, T the gas temperature, e the emissivity and 
p the gas mass density which is related to the electron den- 
sity according to n p — 0.82n e = p/ (2.2\pm p ) (p = 0.6 is 
the mean molecular weight). The equations take into ac- 
count the possible presence of an inflow or outflow and the 
flow velocity v is taken positive outwards. 

The three equations (5) are respectively the mass, mo- 
mentum and energy conservation equations. They have 
been derived (see Sarazin 1988) assuming a steady state; 
the second equation reduces to hydrostatic equilibrium, eq. 
(4), for v — 0. We include in the right-hand-side of the last 
equation, the radiative cooling e, and a heat contribution 
which includes two parts: the thermal conduction e con( i 
which will be widely discussed in the next paragraph, and 
a generic extra-heating term H, which will be studied in 
detail in Section 5.2. 

5.1. Radiative cooling and conduction 

One obvious heating source is the thermal conduction 
which operates when temperature gradients occur, and 
which can have a relatively large efficiency (a fraction of 
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the Spitzer conduction) even in presence of magnetic fields 
(Gruzinov 2002; Narayan and Medvedev 2001). 

Neglecting any extra- heating source (H = in cq. 5), 
and under the assumption of a spherical, steady state, iso- 
baric cooling flow, the last equation of (5) can be rewritten 
in the form: 



M d ( 5kT 



Airr 2 dr \2pm p 



£ &cond ■ 



(6) 



where e is the emissivity, p — 0.6 is the mean molecu- 
lar weight and M is the usual mass deposition rate of the 
cooling flow. 

The heating due to thermal conduction e con( i is given 
by: 



1 d 
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where k is the conductivity. For a highly ionized plasma, 
k is given by the Spitzer (1962) formula: 



K = KS = 



1.84 x 10" 5 (T/°K) 
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(8) 

where In A ~ 40 is the usual Coulomb logarithm. 

Starting from these equations we can derive the conduc- 
tivity k required in M87 to stop the cooling flow (M = 0). 
The inferred values for k are plotted in Fig. A7 (filled cir- 
cles) . We plot k as a function of the temperature in order 
to compare our results with Voigt et al. (2002) and Voigt 
and Fabian (2004) who have performed a similar calcula- 
tion on a set of clusters. We recall that the temperature 
grows with the radius. Hence, the behavior of func- 
tion of the temperature is similar to the behavior of the 
profile of function of the radius. From Fig. A7 we 

can see that the required conductivity has large values for 
small temperatures (i.e. in the central part of the cluster) 
and becomes smaller when the temperature increases, i.e. 
moving towards the outskirts of the cluster. Note that the 
temperature profile in the innermost regions of the clus- 
ter is consistent with being constant. Correspondingly, no 
conduction should be present and the required conductiv- 
ity is consistent with being as high as infinity. Hence, for 
these points, error bars will extent to infinity. In order to 
show in Fig. A7 that the error bar for these points should 
extent to infinity, we plot their error bars with an arrow. 
The solid line in Fig. A7 represents Kg given in eq. (8) and 
dashed line corresponds to 0.3ks which could be the effec- 
tive conductivity in presence of turbulent magnetic fields 
(Gruzinov 2002). Fig. A7 shows that in M87, the ther- 
mal conduction is able to balance radiative cooling only 
in the outer part of the cluster. For r ^ 10 — 20 kpc the 
conductivity required for conduction to balance the cool- 
ing is between 0.3ks and k$. In the inner ^10 — 20 kpc 
in M87 the heating supplied by thermal conduction is not 
enough and an extra-heating, whatever its source might 
be, is needed. The failure of the thermal conduction in 
the core of the cluster is due to the fact that in these re- 
gions the temperature profile flattens (see Fig. A4b) and 
the conduction decreases substantially. At the same time, 
the innermost regions are those where the X-ray emissiv- 



ity is highest and which mostly require a heating source 
to compensate the radiated energy. 

In a recent paper, Voigt et al. (2002) determined the 
conductivity required for the conduction to balance the 
radiation losses for A2199 whose temperature is similar to 
M87 temperature ranging from ~ 2 to ~ 5 keV. For com- 
parison, we plot in Fig. A7 the Voigt et al. (2002) data for 
A2199 derived by modeling the temperature profile with 
two different prescriptions: a power law (triangles) and a 
more complex functional form (squares) which flattens at 
small and large radii (see eq. (6) in Voigt et al. 2002). 
Because of the larger distance of A2199 (z = 0.0309) only 
a couple of points are within the central 10 kpc. Neverthe- 
less, in agreement with our results, they find that for these 
two central bins some extra- heating is needed. The quality 
of the M87 data set allows to highlight the problem of con- 
duction in the core and to analyze it in greater detail than 
for A2199. Clearly, M87 is a good object to test heating 
models. In Fig. A7 we also plot (open circles) the k values 
obtained for M87 by Voigt and Fabian (2004). Their val- 
ues are in reasonable agreement with ours, although their 
analysis procedure is quite different from ours. In Voigt 
and Fabian (2004) only Chandra data are included limit- 
ing the extension of the deprojected region to the inner 
10 kpc and only the IT xspec - mekal model is used in 
fitting spectra extracted from annuli. 

5.2. Heating for M87 and the "effervescent" heating 

model 

In this Section we consider eqs. (5) in their generic form 
in order to determine the extra-heating 7i required to bal- 
ance the radiative cooling in presence of thermal conduc- 
tion for M87. 

Eqs. (5) can be recasted in the form: 
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(9) 



(10) 



This term includes the variation of the energy (per unit 
volume) due to the outflow/inflow and the work (per unit 
volume) done by the system during the outflow/inflow. 
The mass flow rate M is positive for an outflow and neg- 
ative for an inflow. The last equation of (9) provides the 
heating H necessary to balance the radiative cooling e, in 
presence of thermal conduction and steady outflow. The 
deprojected data T, n e (or p), of M87 can be used to solve 
numerically eqs. (9) for M87, deriving M, v and H, once 
we have fixed the fraction / of the Spitzer conductivity 
(eq. 8) and some assumption has been made on M . We 
fix an / = 0.3 efficiency (see Gruzinov 2002) and we set 
the outflow mass rate M = 1.6M /yr. This M value is 
similar to the asymptotic value that RB02 obtained from 
their simulations for the stable final state of the cluster. 
We will discuss further on different values for M and /. 

As far as the velocity is concerned, for the assumed val- 
ues of / and M, v is smaller than a few tens of km/s, 
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for radii larger than ^1 — 2 kpc. Correspondingly, the 
term including the velocity in the momentum conservation 
equation (the second equation of 9) is significantly smaller 
than the total mass being of the order of 10 5 — 1O 7 M 
versus the 10 10 — 10 12 M Q of the total mass, so its contri- 
bution is negligible. Hence, we can state that the cluster is 
almost in hydrostatic equilibrium and the mass estimates 
reported in Fig. A6, where the term related to the outflow 
is neglected, are not affected. In order to have a significant 
contribution from the velocity term and to alter substan- 
tially the hydrostatic equilibrium, M values as large as 
several tens-few hundreds of M Q /yr are required. 

For the considered values of M, the quantity e* in eq. 
(10) is negligible with respect to the emissivity e. Thus, 
the extra-heating TL and the conduction term e con d are 
completely used to balance the radiative cooling. 

For what concerns TL, the heating required in M87 is 
plotted in Fig. A8 (filled circles). As expected, most of 
the heating is required in the central part of the cluster 
(say in the inner ~ 15 — 20 kpc), where conduction is not 
efficient. The heating due to thermal conduction is plot- 
ted in Fig. A8 (dot-dashed line); in the central 10 kpc, 
where the temperature profile becomes flatter, the ther- 
mal conduction drops to zero, apart from the innermost 
bin (~ 1 kpc) where T falls to very small values (see Fig. 
A4b), with a large error bar. The conduction in this bin is 
£cond = 1.46^53 x 10 _24 erg cm~ 3 s _1 , and is in agreement 
with zero within la. 

The heating model developed in Ruszkowski and Bcgel- 
man (2002, RB02 hereafter) and Bcgclman (2001) includes 
a mechanism for heat injection from the central AGN. The 
mechanism has been called "effervescent heating" . The ra- 
dio source is supposed to deposit some buoyant gas bub- 
bles in the ICM, which do not mix and rise through the 
ICM microscopically. The bubbles should expand doing 
work on the surrounding plasma and converting their in- 
ternal energy in heat. The buoyant outflow contribution 
in the energy conservation equation is accounted for in the 
e* term, while TL describes the heat injection due to the 
adiabatic expansion of the bubbles. 

According to the RB02 model, the heating TL is a func- 
tion of the pressure (and its gradient) and can be expressed 
according to: 
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p is the pressure, p is the central pressure, L the time- 
averaged luminosity of the central source, 7& is the adi- 
abatic index of the buoyant gas and ro the inner heating 
cutoff radius. The term 1 — exp{— r/ro) introduces an inner 
cutoff which fixes the scale radius where the bubbles start 
rising buoyantly in the ambient plasma. TL is normalized 
in such a way that, when integrated over the whole cluster, 



the total injected power corresponds to the time-averaged 
power output of the AGN. TL has been derived in eq. (11) 
assuming a steady state for the bubble flux. In order to 
assess if this assumption is reasonable we must compare 
the different timescales involved in the effervescent heat- 
ing mechanism. We can suppose that the AGN is inter- 
mittent (RB02; Reynolds and Begelman 1997) and heats 
the ICM through a succession of outbursts. During each 
outburst, the AGN injects a population of bubbles which 
subsequently rise buoyantly. During the "off" periods, the 
bubbles continue their rise heating the cluster atmosphere. 
If the outbursts follow each other on a timescale which is 
short with respect to the rising timescale of the bubbles, 
then the flux of the bubbles reaches a quasi steady state. 
In fact, the ratio t r i se /ti between the rise timescale t r i se 
and the intermittence timescale U gives the number of pop- 
ulations injected within the ICM within a time t r i se . The 
larger this ratio is, the larger the number of bubble popula- 
tions rising into the cluster atmosphere and the mechanism 
approaches the steady state. The radio galaxies are likely 
to be intermittent on a timescale as short as t% <~ 10 4 — 10 5 
yr (RB02; Reynolds and Begelman 1997). In the next sec- 
tion we will see that the risetime is t r i se ~ 10 8 yr or even 
larger. The value of the ratio trisejU is 10 3 — 10 4 or more; 
therefore the assumption of steady state is tenable and the 
released heating may be treated in a time-averaged sense. 

RB02 also include a convection term in eq. (9) which 
we have neglected. The reason for this choice is twofold. 
First of all, the convection must be limited to the inner- 
most regions of the cluster, in order to allow the presence 
of mctallicity gradients in cluster cores (De Grandi and 
Molendi 2001). Most importantly, a negative gradient for 
the entropy is a necessary condition for the onset of the 
convection. In fact, the condition of instability: 



d_ 

dr 



< 



(14) 



(where 7 is the ratio of the specific heats c p /c v and has 
the value 5/3 for a highly ionized gas) must be fulfilled for 
the convection to operate and it is equivalent to requir- 
ing that V(T/n 2 / 3 ) < 0. This condition is not satisfied in 
M87 where we verified that the entropy is an increasing 
function of the radius, as we show in Fig. A9. 

We compare the values inferred for the extra-heating TL 
term in M87 reported in Fig. A8 (filled circles), with the 
predictions from the RB02 model derived according to eq. 
(11), in order to assess whether, for a reasonable choice 
of the parameters, the heat flux required to balance the 
cooling is compatible with the heat injected by the central 
AGN. 

In order to determine the pressure and the pressure gra- 
dient in eq. (11), we use the analytical expressions (2) and 
(3) for n e and T with the best-fit parameters obtained by 
fitting the deprojected electron density and temperature 
profiles. We use eqs. (11) - (13) as fitting functions for the 
extra-heating, where 7?,, rg and the total normalization A: 



(15) 



are the free parameters. The solid line in Fig. A8 is the 
derived best fit. The model proposed by RB02 seems to 
provide a fair description of the heating needed to balance 
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the cooling flow in the inner <~ 15 — 20 kpc of the clus- 
ter, where conduction is not sufficient. The line seems to 
follow adequately the behavior of the data points. Nev- 
ertheless, the derived fit values have large errors. By 
fixing 7b = 4/3 (which is the adiabatic index for rel- 
ativistic bubbles), we find r = 4.39 ± 1.83 kpc and 
A = (8.35 ± 1.65) x lCT^crgs- 1 (the dashed curve in 
Fig. A8). By using eq. (15), we can also derive the cen- 
tral AGN luminosity L = 5.95 x 10 42 ergs -1 required to 
stop the cooling. 

The three external points in Fig. A8 seem to re- 
quire an additional extra heating, showing an excess with 
respect to the general behavior of the data at large 
radii and with respect to the best fit function shape. 
This excess is related to the flattening of the tempera- 
ture profile in the external regions which dampens con- 
duction. However, it must be noted that the values 
of the heating TL required in these three bins are re- 
spectively: 2.32l°;^ x lCT^crgcm^s- 1 , 5.32j^;j$ x 
10 -28 crgcm -3 s -i and 6.56lg;3 4 x 10 -28 ergcm- 3 s -1 and 
are all in agreement with zero within 2 — 3tr. Hence, the 
evidence for the excess is not particularly strong. 

Our estimated radius ro ^4 — 5 kpc, is comparable 
to the extension of the AGN jet as seen in the Chandra 
image of M87 (e.g. Di Matteo et al. 2003; Young et al. 
2002). The fact that ro is of the same order of magnitude 
of the jet is consistent with a scenario where the efferves- 
cent bubbles are generated through the interaction of the 
radio jet with the cluster atmosphere. Understanding the 
precise nature of such interaction will require considerable 
efforts both on the observational and theoretical side. Our 
aim here is simply to note that our fit docs not rule out 
the possibility that the bubbles are generated through the 
interaction of the radio jet with the cluster atmosphere, 
as would have been the case if, for example, the fitting of 
the effervescent model to M87 had returned an ro value 
10 times larger than the one actually measured. The in- 
ferred luminosity value is similar to the luminosity eval- 
uated for the M87 AGN (Owen, Eilek and Kassim 2000, 
OEK hereafter), - 3 - 4 x 10 42 ergs _1 (see also Di Matteo 
et al. 2003). Slightly different values for the luminosity 
could be inferred with different choices of / and M which 
of course provide different best fit parameters values and 
different related luminosities. When a larger / is consid- 
ered, the higher contribution of the thermal conduction 
reduces the amount of heating needed to balance the ra- 
diative cooling. Setting the conductivity to the Spitzcr 
value (/ = 1) we find best fit values which are similar to 
those inferred for / = 0.3 and provide a slightly lower lu- 
minosity (~ 2 — 3 x 10 42 ergs _1 ). On the contrary, when 
smaller values for / are considered, some additional heat- 
ing is necessary. For / = 0.1 we derive an AGN luminos- 
ity ~l- 2 x 10 43 ergs -1 , which is somewhat larger than 
the OEK estimations. However, for such small efficiencies, 
the shape of the RB02 model no longer provides a good 
description of the data, especially in the central regions. 
Thus, if the contribution of the thermal conduction is too 
small, the "effervescent heating" model is not suitable to 
describe the heating necessary to balance the cooling flow. 

Some variations are found also for different initial val- 
ues of M. We tried 16 and 0.16 M /yr corresponding to 
10 and 0.1 times the original M value we considered. As 



expected, for large values of M the corresponding AGN 
luminosity is significantly enhanced ( ^ /ewlO 43 ergs -1 ) 
since the central AGN must provide a larger quantity of 
energy to the outflowing bubbles. The variations when 
smaller M are considered, are modest, slightly reducing 
the luminosity to ~ 3 — 4 x 10 42 ergs" 1 . 

Note that it is not necessary for the AGN luminosity, ob- 
tained by requiring that the "effervescent heating" model 
balances the cooling flow to be exactly equal to the AGN 
luminosity derived from radio observations. In fact, the 
required luminosity from the model should be regarded as 
a time-averaged power of the AGN, as the AGN dynami- 
cal times are smaller than the radiative cooling flow scale 
times. One should also keep in mind that only a fraction 
of the total power of the AGN is used to quench the cool- 
ing flow and that the luminosity required from the model 
can be significantly smaller than the real AGN luminos- 
ity. From our analysis, we can infer that the values of L 
derived with different choices of / and M are of the same 
order of estimates by OEK and Di Matteo et al. (2003). 

While the luminosity is slightly affected for different 
choices of M and /, ro variations are quite modest and 
the inferred values of the scale radius arc always of the or- 
der of r <~ 4 — 5 kpc, which approximatively corresponds 
to the AGN jet extension. 

5.3. Discussion 

Starting from the results inferred in the previous sec- 
tion, we can try to draw a more general picture, using also 
informations coming from radio observations of M87. 

We can suppose that the buoyant bubbles are radio 
bubbles filled with magnetic field and relativistic parti- 
cles (Gull and Northover 1973; Churazov et al. 2000, 2001; 
Bruggen and Kaiser 2001; Churazov et al. 2002) responsi- 
ble for the synchrotron emission in M87. As already out- 
lined by OEK, the radio structures are highly filamented. 
This suggests that the dimensions of the bubbles are small. 
Enfilin and Heinz (2002) discussed the dynamics of the rise 
of buoyant light bubbles within the cluster atmosphere (see 
also Churazov et al. 2001; Kaiser 2003). The buoyant bub- 
ble rapidly reaches a terminal velocity Vb- In the limit of 
small bubbles, Vb can be estimated by balancing the buoy- 
ancy force with the ram pressure (drag force) of the cluster 
gas. 

The buoyancy force is 

F b = Vg(p - Pb ) = VgpA , (16) 

where V = 4/37rr 3 is the volume of the bubble, n, is the 
bubble radius, g — GM(< r)/r 2 is the local gravity accel- 
eration at the radius r (M is the gravitational mass within 
the radius r); p and pb are respectively the density of the 
ICM and of the bubble. A = (p — pb) / p is the density 
contrast. 

The drag force for subsonic motion can be approximated 

by 

Fdrag = Curlpvl (17) 

with the drag coefficient C ~ 0.5. 

By equating cqs. (16) and (17) the velocity of the bub- 
ble can be determined. Enfilin and Heinz (2002) derived 
Vb under the assumption that the density of the ICM is 
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well described by an isothermal /3-model and the density 
contrast A ~ 1. In this case Vb is a fraction (cx \frbjr~ c ; 
r c is the core radius of the cluster) of the sound velocity. 
Correspondingly, for small bubbles, Vb is subsonic which 
is consistent with the fact that no shocks are detected in 
M87. 

We derived the rise velocity of the bubbles in M87 at a 
radius r ~ 10 kpc, using the results from the deprojection 
for M(< r) and p. The density contrast A can be in- 
ferred considering that the bubbles filled with relativistic 
plasma are in pressure equilibrium with the ICM. Indeed, 
the pressure equilibrium between radio bubbles and ther- 
mal plasma implies that n e kT ~ (3n er (E e ) (where (3 is 
a factor of the order of the unity which accounts for the 
contribution to the pressure of the magnetic field which 
is almost in equipartition with the particle energy). (E e ) 
is the mean energy of the relativistic particles and They IS 
the numerical density of the relativistic particles. The rel- 
ativistic particles mean energy is much larger than that 
of the thermal electrons in the ICM, leading to a density 
contrast A <~ 1. The values of Vb for the bubble radius r& 
varying in [0.01 — 1.] kpc range are plotted in Fig. A10 
and the motion is indeed subsonic. 

It is worth considering that, in a recent analysis of these 
data (Molcndi 2002, M02 hereafter), we have found that, 
cospatially with the radio lobes, there exists a thermal 
component with T ~ \Tjcm- This component is likely 
structured in blobs with typical scales smaller than ~ 100 
pc. The filling factor of these blobs has been estimated to 
be of the order of the percent. 

Some informations about the filling factor of the radio 
bubbles can be obtained using recent results from radio ob- 
servation of M87. Using the standard minimum pressure 
analysis (e.g. Pacholczyk 1970; Burns, Owen and Rudnick 
1979; O'Dea and Owen 1987), OEK evaluate the mag- 
netic field in the lobes and in the filaments visible in the 
radio map and evaluate the minimum pressure of the mag- 
netic field and relativistic particles. They assume that the 
proton-to-electron energy is k — 1 and that £0=1 where 
(f> is the filling factor of the relativistic particles and of the 
magnetic field and £ is the ratio of the true magnetic pres- 
sure (comprehensive of the tension) to the magnetic pres- 
sure B 2 /8tt. For a tangled magnetic field £ = 1/3. The 
estimations for the minimum pressure derived by OEK can 
be compared with the thermal pressure that we can infer 
from our deprojected T and n e profiles. In agreement with 
OEK, we find that the pressure of the relativistic plasma 
is smaller than the thermal pressure by about an order of 
magnitude. By keeping the condition k = 1 and assum- 
ing £ = 1, we can derive the filling factor which reconciles 
the minimum pressure P m i n with the thermal pressure of 
the plasma, considering that the minimum pressure scales 
with (j) according to: 



Pmin{4> = 1) 



= {4>) 



-4/7 



(18) 



The discrepancy of a factor ~ 10 between the two pres- 
sures can be eliminated by assuming a filling factor ~ 1%. 
By using the estimations of P m in from OEK in different 
plasma of the radio lobes and filaments, we infer filling fac- 
tors of few %. Assuming £ = 1/3 (tangled magnetic field) 
reduces by a factor ~ 2. The filling factors of the ra- 



dio bubbles and of the cold thermal blobs are of the same 
order. 

The survival of the cold thermal blobs in the hotter 
ICM requires that thermal conduction be substantially 
suppressed; this may happen if these blobs are tied to 
the radio bubbles and magnetic fields shield them from 
collisions with ICM particles. Since the filling factors of 
blobs and bubbles are similar and the density of the blobs 
is about twice that of the ICM (M02), the mean density of 
each bubble+blob is about that of the surrounding ICM. 
Thus, assuming that the blobs are tied to the bubbles and 
that they occupy similar volumes, their density contrast 
with the surrounding ICM will be small. In Fig. A10 the 
dotted line plots the bubble rise velocity Vb for a A = 0.1. 
Starting from the rise velocity, it is straightforward to 
derive the rise timescale of the bubbles. For a density 
A <~ 1 and a bubble radius of <~ 100 pc, the risetime is 
t-rise ~ 10 8 yr at r ~ 10 kpc and it becomes even larger 
for smaller bubbles and smaller density contrasts, holding 
t rise oc (r&A) ' . As outlined previously, this value is sig- 
nificantly larger than the duty cycle timescale ti of inter- 
mittency of the AGN, so that the ratio t r i se /ti £ 10 3 — 10 4 
and the mechanism approaches the steady state. 

It is worth noting that the above picture refers to radio 
bubbles and cool X-ray blobs located in the lobes. How- 
ever, the heating process related to the rise of the bubbles 
must be isotropic throughout the cluster in order to bal- 
ance the cooling flow. Under the assumption that the radio 
galaxy is intermittent on a short timescale, the mechanism 
is expected to heat isotropically the ICM since it is likely 
that no direction is preferred for the AGN ejection; we also 
recall that when the AGN is turned off the bubbles con- 
tinue their rise within the ICM heating it up; within this 
picture, by averaging on a cooling time, the radio bubble 
populations arc likely to be isotropically distributed in the 
cluster. 

Nevertheless, the cool X-ray blobs are detected only in 
the lobes regions. So the picture emerging here is that 
of an AGN which injects radio bubbles, in all the direc- 
tions, through a succession of outbursts. In the lobes, the 
outburst is occurring at the present time and, here, also 
the cool X-ray blobs are present. They are tied to the 
radio bubbles and the magnetic fields shield them from 
collisions with ICM particles. The thermal conduction is 
inhibited allowing the cool blobs to survive in a hotter 
ambient medium. In the halo, where the radio bubbles 
have been injected during a past activity of the AGN, only 
"old" populations of bubbles which are buoyantly rising 
are present. It is likely that the cool blobs which were tied 
to the radio bubbles when the AGN was active there, have 
thermalized, the magnetic fields slow down but do not stop 
entirely thermalization so that only a single phase gas is 
detected in the halo regions. 

6. SUMMARY 

The crisis of the standard cooling flow model brought 
about by Chandra and XMM-Newton observations of 
galaxy clusters, has led to the development of several heat- 
ing models with the aim of identifying a mechanism able 
to quench the cooling flow. 

We have used observations of Virgo/M87 to inspect 
the dynamics of the gas in the center of the cluster, 
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and to study the heating processes able to balance radia- 
tive losses. We combined the observations of three satel- 
lites, namely XMM-Newton, Chandra and Beppo-SAX. By 
means of the spectral deprojection technique, we inferred 
the profiles for the temperature T and the electron density 
n e for M87. The temperature profile drops by a factor of 
<~ 2 in the inner ( ,$ 15 kpc). The electron density n e is 
well described by a 2-/3 profile. Starting from these pro- 
files, we derived some related physical quantities such as 
the gravitational mass and the entropy profile. The grav- 
itational mass profile shows a plateaux at <~ 10 — 15 kpc, 
which could correspond to the edge of the central cD. 

In agreement with Voigt et al. (2002) and Voigt and 
Fabian (2004) results on a set of cluster, we found that 
the thermal conduction in M87 can balance the cooling 
flow only in the outer part of the core (say, r £ 15 kpc), 
while in the inner 10 — 15 kpc, where the temperature 
profile becomes flatter and the conduction is no longer 
efficient, some extra-heating is required. We have deter- 
mined the extra-heating needed to balance the cooling flow 
in M87, assuming a thermal conduction efficiency / = 0.3 
with respect to the Spitzer value, and an outflow mass rate 
M = 1.6M©/yr. The high quality of our combined dataset 
allows us to inspect properly the innermost regions, deriv- 
ing an accurate profile for the extra-heating term H in 
the regions where thermal conduction is not sufficient to 
quench the cooling flow. 

Several models and simulations concerning heating 
mechanisms through buoyant gas in the cluster ICM have 
been recently proposed in the literature (e.g. Ruszkowski 
and Begelman 2002; Begclman 2001; Churazov et al. 
2002; Bohringer et al. 2002; Briiggen and Kaiser 2002a,b; 
Briiggen et al. 2002). We have assumed that the heat- 
ing is provided by the central AGN by means of depo- 
sition of buoyant bubbles in the ICM according to the 
model proposed by Ruszkowski and Begelman (2002). The 
bubbles rise through the ICM and expand doing work 
on the surrounding plasma and heating it up. We fit- 
ted the extra-heating required in M87, with the heating 
functions proposed by RB02 (see eqs. 11 - 13). The 
RB02 model seems suitable to describe the behavior of 
the data. By fixing the adiabatic index 75 = 4/3 (relativis- 
ts bubbles), we find a scale radius r = 4.39 ± 1.83 kpc 
and A = (8.35 ± 1.65) x 10 _24 ergs _1 , which corresponds 
to a central AGN luminosity L = 5.95 x 10 42 ergs _1 . 
The scale radius is of the order of the extension of the 
AGN jet and the inferred AGN luminosity is similar to 
the one estimated by Owen, Eilek and Kassim (2000) 
and Di Matteo et al. (2003) L ~ 3 - 4 x 10 42 erg/s _1 . 
Smaller conductivity efficiencies (/ = 0.1) or larger out- 
flow mass rates (M = 16M Q /yr) provide AGN luminosi- 
ties (L <~ fewlO 43 erg s _1 ) which are somewhat higher 
than the estimations of Owen, Eilek and Kassim (2000) 
and Di Matteo et al. (2003). However, if the efficiency of 
the thermal conduction is reduced (/ = 0.1) the model 
functions of RB02 seem no longer suitable to describe the 
heating needed to balance the cooling flow. For higher con- 



duction efficiencies (/ = 1) or smaller outflow rates values 
(M = 0.16M Q /yr) the inferred luminosity is slightly re- 
duced (L ~ 3 — 4 x 10 42 erg/s -1 ). The different values 
derived for L arc of the same order of the luminosity mea- 
sures for the AGN luminosity in M87 obtained by Owen, 
Eilek and Kassim (2000) and Di Matteo et al. (2003). In 
all the cases considered, the inferred scale radius r , which 
fixes the radius at which bubbles are deposited and start 
rising, is ro ~ 4—5 kpc which approximativcly corresponds 
to the AGN jet extension. 

Finally, we discussed a scenario where the bubbles are 
filled with rclativistic particles and magnetic field, respon- 
sible for the radio emission in M87. In this scenario the 
density contrast A is expected to be as large as 1. The 
buoyant velocity Vf, of the bubbles can be derived by bal- 
ancing the buoyant force with the drag force. For small 
bubbles, the rise velocity is subsonic. Under the hypothesis 
of equipartition between relativistic particles and magnetic 
field in the bubbles and of equilibrium pressure with the 
thermal ICM, we evaluated the filling factor (f> of the radio 
bubbles. We find <f> ~ 0.01 which is of the same order of 
the filling factor of the cool thermal component observed 
in the regions of the radio lobes (M02). Hence, we suggest 
that this cool thermal component is structured in blobs 
tied to the radio bubbles. The thermal conduction, which 
should rapidly thcrmalizcd the cool blobs, is suppressed by 
the magnetic fields of the radio bubbles. The density con- 
trast of the buoyant bubble+blob system is A < 1 further 
reducing the rising velocity. 

The radio galaxies are likely to be intermittent on a 
timescale t{ ~ 10 4 — 10 5 yr and they are supposed to 
heat the ICM through a succession of outbursts. Dur- 
ing each outburst a population of radio bubbles is injected 
into the ICM. The bubbles rise buoyantly in the intraclus- 
ter gas, heating it up. The outbursts follow each other 
on small timescale t{ which is much shorter than the rise 
timescale t rise £ 10 8 yr of the bubbles so that the mech- 
anism is isotropic throughout the cluster and approaches 
the steady state. 

The X-ray cool blobs are detected in the radio lobes 
where the injection of radio bubbles is occurring at the 
present time. The blobs are tied to the radio bubbles so 
that the thermal conduction is highly suppressed by the 
magnetic fields. In the radio halo, the radio bubbles injec- 
tion occurred in the past. The radio bubbles are buoyantly 
rising in the cluster atmosphere and it is likely that the 
X-ray cool blobs which were tied to the bubbles when the 
AGN was active in those directions, have in the meanwhile 
thcrmalizcd so that only a single phase thermal component 
is present here. 
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APPENDIX 



CORRECTION FACTOR IN THE DEPROJECTION RECIPE FOR THE BACKGROUND EMISSION 

By using eq. (1) the deprojected variable f shell is obtained from F ring by subtracting off the contribution of the outer 
shells, starting from the outermost annulus and moving inwards. This basic prescription assumes that there is no emission 
beyond the maximum radius R m ax- Thus, the basic recipe (1) must be corrected to account for a contribution to the 
X-ray spectra from the gas beyond R ma x- 

In practice, in eq. (1) the 2-D variable to be deprojected F ring is replaced by an effective one F^fJ g defined by: 

^ring ~ ^ring 9ring ' F n • A r i ng /A n , (-^-1) 

where A ring is the area of the ring, F ring the variable to be deprojected and g r i ng the correction factor of the ring. A n 
and F n are the area and the variable to be deprojected of the outer ring. 

In order to determine g r ing, some assumptions on the shape at large radii for the different deprojected quantities must 
be made. The classic solution consists in assuming that all eT s h e ii, £ and n e have an / oc r~ a dependence at large radii, 
with a = 4. For this standard assumption, the correction factor gj can be determined analytically (McLaughlin 1999). 

Hence, the correction to the j-th ring for the contribution coming from the outer part of the cluster is: 



Rl-Rl_J R L dbb $V 



+°° dz 



Rl-b 2 (b 2 +z 2 ) 



90 ini^ir-^ ' 



(b2 +z 2)a/2 



R and b refer to the 2-D radii of the rings and z is the line-of-sight integration variable and / oc r~ a has been used. Of 
course r 2 = b 2 + z 2 holds. 

However, the standard assumption a = 4 is quite simplistic and holds only for an isothermal cluster with n e oc r~ 2 at 
large radii (e.g. in the usual (3- model with [3 = 2/3). On the contrary, we prefer considering the dependence r~ a with a 
generic a, different for each deprojected quantity. While for a = 4 the correction factor could be calculated analytically, 
for a generic a it must be determined numerically. We evaluated the gj factor of eq. (A2) by truncating the integrals at 
a 10R max external radius, with steps 0.01R max wide. We verified that the numerical method for the simple case a = 4 
provides results which differ by no more than a few percent from the analytical values, and in any case the effect is always 
limited to the outermost bins. In order to find out the correct a for each quantity to be deprojected, we applied itcratively 
the deprojection. We started from reasonable values of a. Then we applied the deprojection and derive the new a values 
from the slopes of the deprojected profiles at large radii. We used these new a values to determine the correction factor 
of eq. (A2) and applied again the deprojection working out new a values. We stopped when all the new a values differed 
from the starting values by less than 4%. 
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Fig. Al. — Emission weighted temperatures for the single phase gas of the IT regions (open circles) and for the hot component (triangles) 
of the 2T regions. Error bars are plotted only for few representative points. These points are plotted at slightly larger radii for a clearer view 
of the error bars amplitude. 
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Fig. A2. — (a) Emission-weighted temperature profile for M87 obtained from the spectral analysis described in §2; (b) Normalized Emission 
Integral (NEI) profile per unit area; NEI is given in XSPEC units, i.e. NEI = ,}° , vi EI, where d ang is the angular distance of M87 in 

ang\ ' z l 

cm, z the redshift and El in cm . The Area is in arcmin . In both panels, circles refer to XMM-Newton data, squares to Chandra data 
and triangles refer to Beppo-SAX data. 
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Fig. A3. — Emission-weighted temperature profile for M87 obtained from the spectral analysis described in §2, with the Chandra and 
Beppo-SAX data sets renormalized in order to match the XMM-Newton profile. Symbols are the same as in Fig. A2. 
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Fig. A4. — (a) Electron density derived with the deprojection method (filled circles). The open diamonds represent the n e profile after the 
smoothing operation. The solid line is the best fit profile where a 2-/3 model (see eq. 2) has been used, (b) Deprojected temperature profile 
versus radius (filled circles). The open diamonds represent the T profile after the smoothing operation. The solid line is the best fit where 
the expression given in eq. (3) has been used. Error bars in both the panels have been obtained by 1000 Monte Carlo simulations on initial 
EI and Tew ■ The details on the smoothing operation and on the determination of the error bars are discussed in §4.1. The triangles are the 
deprojected density and temperature profiles derived by Matsushita ct al. (2002). 
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Fig. A5. — Deprojected density profile (a) and deprojected temperature profile (b) versus the radius for the whole cluster (filled dots) and 
for the non-radio regions (open diamonds). No significant differences are evident. 
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Fig. A6. — Gravitational mass derived through the hydrostatic equilibrium equation. The filled circles refer to the profile inferred 
by smoothing the temperature and the density profiles. The open diamonds plot the gravitational mass derived without any smoothing 
operation on n e and T, errors have been derived using the standard Monte Carlo technique. Three of these points have values near to zero 
and only the upper limit of their error bar has been plotted here. The solid curve is the analytical mass obtained using the best fit profiles 
for n e and T (eqs. 2 and 3). For comparison, the grey-shaded regions report the gravitational mass derived by Nulsen and Bohringcr (1995) 
using ROSAT-PSPC data. 
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Fig. A7. — The conductivity coefficient required for conduction to balance radiation losses, for M87 (close circles). The points which have 
an error bar which extents to infinity have been plotted with an arrow. The solid line is the Spitzer conductivity and the dashed line is one 
third of the Spitzer conductivity. The open circles refer to results for M87 derived from Voigt and Fabian (2004). The triangles and the 
squares refer to Voigt et al. (2002) data for A2199 derived with two different prescriptions (see the text for details). 
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Fig. A8. — The heating required (filled circles) to balance radiation losses in M87. The dot-dashed line is the heating due to thermal 
conduction. The solid line is the best fit obtained fitting the data set with the RB02 model. The dashed line is the best fit obtained fixing 
76 = 4/3. 



22 



Ghizzardi et al. 



2.0 



> 



CjjO 
O 



1.5 - 



1.0 - 



0.5 



it* 



10 

r (kpc) 



100 



Fig. A9.— Entropy profile for M87. 




. A10. — Rise velocity for the buoyant bubbles for density contrast between the bubble and the ambient ICM A 
0.1 (dashed line). All the quantities have been evaluated at r ~ 10 kpc. 
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Table Al 

Emission-weighted temperatures in keV and Normalized Emission Integral (NEI) per unit area; NEI is given in 
xspec units, i.e. NEI = - J° . yi EI, where dang is the angular distance of M87 in cm, z the redshift and EI in 

cm -3 . The Area is in arcmin 2 . 
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Table A 2 

Temperature, electron density and emissivity values obtained with the deprojection technique for M87. 
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